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Abstract 

In  this  contribution,  we  report  on  some  results  recently  obtained  in 
joint  work  with  Maksymilian  Dryja.  We  first  study  an  additive  variant 
of  Schwarz'  alternating  algorithm  and  establish  that  a  fast  method  of 
this  kind  can  be  devised  which  is  optimal  in  the  sense  that  the  number 
of  conjugate  gradient  iterations  required,  to  reach  a  certain  tolerance, 
is  independent  of  the  mesh  size  as  well  as  the  number  of  subregions.  An 
interesting  feature  of  this  algorithm  is  that  all  the  subproblems  can  be 
solved  at  the  same  time.  The  method  is  therefore  quite  promising  for 
parallel  computers.  Using  a  similar  mathematical  framework,  we  then 
consider  the  solution  of  elliptic  finite  element  problems  on  composite 
meshes.  Such  problems  can  be  built  up  systematically  by  introducing 
a  basic  finite  element  approximation  on  the  entire  region  and  then  re- 
peatedly selecting  subregions,  and  subregions  of  subregions,  where  the 
finite  element  model  is  further  refined.  We  consider  conjugate  gradient 
algorithms  where,  in  each  iteration,  problems  on  the  subregions  rep- 
resenting finite  element  models  prior  to  further  refinement  are  solved. 
This  makes  it  possible  to  use  solvers  for  problems  with  uniform  or  rel- 
atively uniform  mesh  sizes,  while  the  composite  mesh  can  be  strongly 
graded.  We  remark  that  this  work  is  technically  quite  closely  related 
to  our  previous  work  on  iterative  substructuring  methods,  which  are 
domain  decomposition  algorithms  using  non-overlapping  subregions. 


1     Introduction 

In  this  paper,  we  report  on  results  obtained  in  joint  work  with  Maksymil- 
ian  Dryja  of  the  University  of  Warsaw,  Poland.  Our  work  began  with  the 
study  of  a  recent  paper  by  P.-L.  Lions  in  which  a  variational  framework 
is  introduced  for  the  classical,  multiplicative  (sequential)  Schwarz'  method; 
see  Schwarz  [19]  and  Lions  [13].  As  pointed  out  by  Lions,  the  variational 
framework  is  far  from  new;  cf.  [20],  [2],  [1]  and  other  references  given  in 
Lions  [13]. 

In  the  second  section  of  this  paper,  we  introduce  a  framework  similar 
to  Lions'  and  also  an  additive  (  parallel)  version  of  the  Schwarz  algorithm. 
While  Lions  considered  the  continuous  problems,  we  work  consistently  with 
conforming  finite  element  approximations;  cf.  [6]. 

In  view  of  the  interest  in  parallel  computers  with  many  processors,  we 
are  principally  interested  in  the  case  of  many  subregions.  Our  main  result 
is  that  we  can  design  additive  methods  of  this  kind  which  have  rates  of 
convergence  which  are  independent  of  the  number  of  subregions  and  the 
number  of  unknowns.  No  proofs  are  given  in  this  paper;  cf.  Dryja  and 
Widlund  [8]  for  details.  We  note  that  the  technical  tools  used  there  are 
similar  to  those  of  our  previous  work  on  so-called  iterative  substructuring 
methods,  which  are  domain  decomposition  methods  with  nonoverlapping 
subregions;  cf.  e.g.  [7]  and  [21]  and  the  articles  cited  in  those  papers. 

In  the  final  section,  we  discuss  optimal  iterative  refinement  methods. 
This  work  begun  with  the  discovery  that  the  so  called  FAC  and  AFAC 
methods,  which  have  been  studied  by  McCormick  and  his  coworkers  [15], 
[14],  [16],  [17]  and  [12],  and  Bramble,  Ewing  et  al.  [5],  [9],  are  structured 
quite  similarly  to  the  classical  Schwarz  procedure  and  its  additive  variant, 
respectively.  In  our  view,  the  central  theoretical  issue  of  the  iterative  re- 
finement algorithms  is  the  design  and  study  of  algorithms  for  which  the 
rate  of  convergence  is  independent  of  the  number  of  subproblems,  i.e.  the 
number  of  refinement  levels,  as  well  as  the  mesh  size.  After  introducing  cer- 
tain finite  element  models  on  composite  meshes,  which  can  systematically 
be  built  up  inside  a  framework  of  conforming  finite  elements,  we  describe 
several  algorithms,  for  which  we  have  established  optimality  in  this  sense; 
cf.  [23]. 

Let  us  also  note  that  we  have  discovered  that  the  very  interesting  so 
called  hierarchical  basis  method,  introduced  recently  by  Yserentant  [24], 
can  be  derived  and  understood  inside  our  general  framework.  We  plan  to 
return  to  this  topic  elsewhere. 


Since  the  theme  of  this  seminar  is  boundary  integral  equations  and 
boundary  element  methods,  a  remark  on  a  connection  between  this  main 
topic  and  domain  decomposition  algorithms  might  be  in  order.  In  our  early 
work  on  domain  decomposition  algorithms  for  regions  divided  into  two  non- 
overlapping  subdomains,  we  established  that  certain  discrete  boundary  op- 
erators play  a  crucial  role  in  the  description  and  analysis  of  the  so  called 
Neumann-Dirichlet  algorithm;  cf.  Bj0rstad  and  Widlund  [4].  In  the  basic 
form  of  this  algorithm,  the  problems  on  the  different  subregions  are  solved 
one  after  the  other,  using  natural  and  essential  boundary  conditions  in  al- 
ternate steps.  The  resulting  iteration  matrix  can  be  represented  in  terms  of 
so  called  Schur  complements  of  the  stiffness  matrices  of  the  subregions.  The 
Schur  complements  are  obtained  by  eliminating  all  the  interior  variables  in 
the  subregions.  We  can  equally  well  regard  them  as  the  restrictions  of  the 
quadratic  forms,  which  represent  the  energy  of  the  subregions,  to  6pa<;es 
of  discrete  harmonic  functions.  These  boundary  operators  map  the  trace 
of  the  finite  element  functions  onto  a  natural  approximation  of  the  normal 
derivative  at  the  boundary  of  the  corresponding  discrete  harmonic  function. 
By  using  finite  element  techniques,  we  have  established  a  bound,  which  is 
a  proper  analog  of  a  standard  estimate  of  elliptic  theory.  This  bound  is 
crucial  in  the  proof  that  the  rate  of  convergence  of  the  iterative  scheme  is 
independent  of  the  mesh  size. 

These  operators  can  be  regarded  as  discretizations  of  what  most  appro- 
priately have  been  called  the  canonical  integral  operators;  cf.  Feng  [10]. 
The  same  bound  is  available  for  these  boundary  integral  operators  and  their 
discretizations  as  for  the  finite  element  operator.  We  can  therefore  apply 
the  Neumann-  Dirichlet  algorithm  to  coupled  finite  element  -  boundary  in- 
tegral equation  models.  In  particular,  if  the  same  boundary  space  is  used 
for  the  two  subproblems,  we  can  establish  that  the  rate  of  convergence  of 
the  resulting  method  is  independent  of  the  number  of  degrees  of  freedom. 


2     An  Additive  Variant  of  the  Schwarz  Alternat- 
ing Method 

We  consider  linear,  self  adjoint,  elliptic  problems  discretized  by  finite  element 
methods  on  bounded  Lipschitz  regions.  We  assume  that  fi  is  a  region  in  J?", 
that  the  differential  operator  is  the  Laplacian  and  that  we  use  continuous, 
piecewise  linear  finite  elements.  The  theory  can  equally  well  be  developed 
for  any  linear  elliptic  problem  which  can  be  formulated  as  a  minimization 
problem,  and  we  can  also  use  arbitrary  conforming  finite  elements.  The 
continuous  and  discrete  problems  are  of  the  form 

a{u,v)  =  f{v),W  V  e   V, 
and 

a{UK,VH)  =  f(VH),'^VHe    V\  (1) 

respectively.  The  bilinear  form  is  defined  by 


a(u,v)  =   /  Vu-Vvdx 


It  defines  a  semi-norm  \u\ffi  =  {a{u,u)y'^  in  the  Sobolev  space  n^{fl  ). 
For  simplicity,  we  assume  that  we  have  zero  Dirichlet  boundary  conditions. 
All  elements  of  V  and  its  subspace  V^  vanish  on  5fi,  the  boundary  of  fi  . 
We  introduce  the  triangulation  of  Cl  in  the  following  way.  We  first  divide  the 
region  into  non-overlapping  substructures  fii ,  t  =  1,  •  •  • ,  JV  .  To  simplify  the 
description  of  our  algorithm,  we  confine  our  study  to  triangular  (simplicial) 
substructures.  In  such  a  case,  the  original  region  must  of  course  be  a  polygon 
(polyhedron).  We  also  adopt  the  common  assumption  in  finite  element 
theory  that  all  substructures  are  shape  regular  in  the  sense  that  Hi,  the 
diameter  of  the  fi,  divided  by  the  radius  of  the  largest  inscribed  sphere 
in  fli  is  uniformly  bounded.  AU  the  substructures  fi,  are  divided  into 
elements.  The  elements  are  also  shape  regular  in  the  same  sense.  Since 
Schwarz-type  domain  decomposition  algorithms  use  overlapping  subregions, 
we  extend  each  substructure  to  a  larger  region  il[.  We  assume  that  the 
distance  between  the  boundaries  5Q,  and  dil;  is  bounded  from  below  by 
a  fixed  fraction  of  Hi,  and  that  dfl\  does  not  cut  through  any  element.  We 
make  the  same  construction  for  the  substructures  that  meet  the  boundary 
except  that  we  cut  off  the  part  of  Ct'-  that  is  outside  of  ft. 

We  remark  that  some  of  the  more  intricate  issues  in  the  analysis  of  the 
Schwarz  methods  arise  when  the  boundaries  of  the  different  subdomains  QJ 


intersect  at  one  or  several  points;  cf.  the  discussion  in  Lions  [13].  This 
happens  for  example  if  the  region  is  L-shaped  and  is  partitioned  into  two 
overlapping  rectangles.  We  discuss  this  matter  in  our  more  complete  paper 
on  this  subject. 

Our  finite  element  space  is  represented  as  the  sum  of  N+1  subspaces 

V''  =  Vo''  +  Vl"  +  ...  +  V;^ . 

The  first  subspace  Vq''  which  we  also  caU  V^,  is  special.  It  is  the  space 
of  continuous,  piecewise  linear  functions  on  the  coarse  mesh  defined  by  the 
substructures  Cli.  We  use  this  spa<;e  to  provide  a  mechanism  of  global  trans- 
portation of  information  similar  to  the  role  of  coarse  grid  solvers  in  miiltigrid 
methods;  cf.  [8],  [21]  or  [22].  It  is  quite  easy  to  show  that  without  such  a 
device,  any  additive  Schwarz  method  would  converge  no  faster  than  the  con- 
jugate gradient  method  applied  directly  to  the  same  coarse  finite  element 
discretization.  The  other  subspaces  are  related  to  the  subdomains,  in  the 
same  way  as  the  traditional  Schwafz  algorithm. 

For  self-adjoint  elliptic  problems,  the  finite  element  solution  is  the  pro- 
jection of  the  exact  solution  onto  a  finite-dimensional  space.  Different  vari- 
ants of  Schwarz'  algorithm,  in  particular  the  additive  form  considered  in 
this  paper,  can  also  conveniently  be  described  in  terms  of  projections  onto 
subspaces  related  to  overlapping  subregions  which  cover  the  region.  These 
projections  are  orthogonal  with  respect  to  the  symmetric  bilinear  form  as- 
sociated with  the  elliptic  operator.  The  computation  of  the  projection  of  an 
arbitrary  function  onto  V^  involves  the  solution  of  a  standard  finite  element 
linear  system  of  algebraic  equations  which  has  on  the  order  of  N  unknowns. 
This  coarse,  global  approximation  of  the  elliptic  equation  is  of  the  same  type 
as  the  local  problems  associated  with  the  subdomains.  It  is  often  feasible 
to  handle  all  these  linear  systems  of  equations  by  a  direct  method  such  as 
Gaussian  elimination.  The  only  real  difference  between  the  problem  related 
to  the  first  subspace  and  the  others  lies  in  the  way  the  right-hand  side  of  the 
linear  system  is  generated  as  weighted  averages  with  weights  determined  by 
the  basis  functions  associted  with  the  coarse  mesh.  We  note  that  if  we  make 
the  dimension  of  all  the  subspaces  approximately  equal,  then  we  have  A'^  -|- 1 
linear  systems,  each  with  about  than  N  unknowns,  to  solve  in  each  step  of 
the  iterative  solution  of  a  linear  system  with  about  N^  unknowns.  In  this 
algorithm  two  mesh  sizes  are  used,  but  we  could  also  introduce  additional, 
intermediary  finite  element  models. 

The  subspaces,  V/^ ,  which  are  related  to  interior  substructures,  are  de- 
fined as  V'*  n  nl{Q'i)  .  The  space  H^    is,  as  usual,  the  subspace  of  H^ 


functions  with  zero  trace.  It  is  well  known  that  this  space  can  be  extended 
continuously  by  zero  to  the  rest  of  Q.  After  such  an  extension,  we  can  regard 
Vj''  as  a  subspace  of  V^  .  We  construct  subspaces  related  to  the  boundary 
substructures  similarly. 

The  projection  Pyh  =  Pi  is  defined,  for  all  of  V^,  by  the  unique  element 
of  V/^,  which  satisfies 

a{PiVh,<i>h)  =  a{vh,<t>h),  ^4>h  e  V>.  (2) 

Lions  [13]  has  shown  that  the  error  propagation  operator  of  the  standard 
multiplicative  variant  of  the  Schwarz  method  can  be  written  as 

(/-P2)(/-Px), 

for  the  case  of  two  subregions.  We  can  therefore  view  that  algorithm  as  a 
simple  iterative  method  for  solving 

(A  +  P2  -  P2Pl)Uh  =  9h, 

with  an  appropriate  right-hand  side  g^  .  We  note  that  the  operator  is  a 
polynomial  of  degree  two  and  thus  is  not  ideal  for  parallel  computing  since 
two  sequential  steps  are  involved.  If  we  use  more  than  two  subspaces  this 
effect  is  further  pronounced,  even  if  the  degree  of  the  polynomial  representing 
the  multiplicative  algorithm  often  is  lower  than  the  maximal.  This  is  so 
because  products  of  projections  associated  with  subregions,  which  do  not 
overlap,  vanish. 

The  basic  idea  behind  the  additive  form  of  the  algorithm  is  to  work 
with  the  simplest  possible  polynomial  in  the  projections.  We  thus  solve  the 
equation 

Puh  =  (Po  +  Pi  +  ....  +  PnW  =  9h  ,  (3) 

by  an  iterative  method.  Since  the  operator  P  can  be  shown  to  be  symmetric 
and  positive  definite,  the  method  of  choice  is  the  conjugate  gradient  method. 
We  must  make  sure  that  equation  3  has  the  same  solution  as  equation  1,  i.e. 
we  must  find  the  correct  right-hand  side.  Since  by  equation  1,  we  have 

a(uh,4>h)  =  fi4>h), 

we  can  construct  the  right-hand  side  g'f^  by  solving  equation  2  for  all  values 
of  i  and  by  adding  the  resulting  vectors.  It  is  similarly  possible  to  apply 
the  operator  P  of  equation  3  to  any  given  element  of  V''  by  applying  each 


projection  P,  once  and  adding  the  results.  Most  of  the  work,  in  particular 
that  which  involves  the  individual  projections,  can  be  carried  out  in  parallel. 
It  is  well  known  that  the  number  of  steps  required  to  decrease  an  appro- 
priate norm  of  the  error  of  a  conjugate  gradient  iteration  by  a  fixed  factor  is 
proportional  to  y/n  ,  where  k  is  the  condition  number  of  P  ;  see  e.g.  Golub 
and  Van  Loan  [11].  We  therefore  need  to  establish  that  the  operator  P 
of  equation  3  is  not  only  invertible  but  that  satisfactory  upper  and  lower 
bounds  on  its  eigenvalues  can  be  obtained.  We  have  proven  the  following 
theorem;  cf.  [8]. 

Theorem  1  k(P)   is  hounded  by  a  constant. 

Here,  and  elsewhere  in  this  paper,  the  constants  are  independent  of  the 
diameters  of  of  the  substructures  and  the  individual  elements  but  they  vary 
with  the  shape  of  the  elements  and  substructures.  Thus,  in  two  dimensions, 
they  depend  on  the  minimum  angle.  As  becomes  apparent  by  examing  the 
proof,  the  constant  in  the  theorem  grows  if  the  overlap  of  the  subregions 
shrinks. 

We  have  discovered  the  additive  algorithm  independently,  but  its  origin 
is  not  clear.  We  are  quite  interested  in  finding  out  about  any  previous  work 
of  this  nature.  We  have  found  an  algorithm  quite  similar  to  ours,  in  a  paper 
by  Nepomnyashchikh  [18]  on  iterative  substructuring  methods  on  regions 
subdivided  into  a  fixed  number  of  subregions. 


3     Iterative  Refinement  Algorithms 

Finite  element  models  on  composite  meshes  can  systematically  be  built  up 
inside  a  conforming  finite  element  framework.  Conforming  elements  allow 
us  to  use  a  number  of  technical  tools  otherwise  not  available.  We  build 
our  models  systematically  by  first  introducing  a  ba^ic  finite  element  ap- 
proximation on  the  entire  region  and  then  repeatedly  selecting  subregions, 
and  subre^ons  of  subregions,  where  the  finite  element  model  is  further  re- 
fined. We  consider  conjugate  gradient  algorithms  where,  in  each  iteration, 
problems  representing  relatively  simple  finite  element  models  on  the  origi- 
nal regions,  and  the  subregions,  prior  to  further  refinement  are  solved.  We 
call  these  the  standard  problems.  We  thus  use  solvers  for  problems  with 
uniform  or  relatively  uniform  mesh  sizes,  while  the  composite  mesh  can  be 
strongly  graded.  As  shown  in  [23],  it  is  also  possible  to  solve  the  standard 
problems  inexcictly.  This  type  of  algorithms  is  of  interest  when  we  wish  to 
use  adaptive  meshes  to  increase  accuracy  locally.  It  is  also  important  to  note 
that  existing  industrial  codes  can  be  improved  relatively  easily  by  adding 
refinement  regions  in  this  manner;  of  [5],  [9]. 

We  study  two  families  of  algorithms  which  are  multiplicative  and  addi- 
tive respectively  in  the  same  sense  as  before.  The  additive  algorithms  are 
particularly  well  suited  to  parallel  computing  since,  in  each  iteration  step, 
we  can  simultaneously  solve  all  the  standard  problems  corresponding  to  all 
the  different  levels  of  refinement.  Synchronization  between  the  processors  is 
required  once  in  each  conjugate  gradient  step  when  the  residual  is  computed 
and  assembled.  One  or  a  group  of  processors  can  therefore  be  assigned  in  a 
straightforward  way  to  each  of  the  standard  subproblems. 

We  introduce  some  notations.  Denote  the  original  region  il  by  Qi  and 
let  fii  C  fit-i ,  i  =  1  ■  •  •  fc,  be  subregions  of  subregions  of  the  original  region. 
As  explained  further  in  [23],  we  impose  certain  mild  geometric  conditions 
on  fl,-i  \  n,.  Let  V^' ,  i  =  1  •  ■  ■  fc,  be  a  family  of  conforming  finite  element 
spaces  such  that  V''-'  fl  HHVli)  C  V^'  C  Hl{Qi).  Our  goal  is  to  solve 
equation  1  using 

The  fundamental  building  blocks  of  our  algorithms  are  the  Pj,i  <  j,  the 
projections  onto  the  spaces  V^'  n^o(^j)-  ^^  note  that  if  j  >  i,  then  we 
solve  a  problem  on  fij  with  a  coarser  mesh  than  if  V^'  were  used.  Such  a 
problem  is  therefore  relatively  less  expensive. 

We  can  now  describe  a  basic  multiplicative  refinement  algorithm.  The 
operator  corresponding  to  this  method  is  symmetric  and  we  can  accelerate 


it  with  the  standard  conjugate  gradient  method.  In  each  of  2fc  —  1  substeps 
of  a  basic  iteration  step,  we  compute  the  residual,  using  the  just  updated 
solution.  The  bilinear  form  corresponding  to  this  residual  now  serves  as 
a  right-hand  side.  By  solving  in  the  subspace  V'*',  we  compute  the  next 
correction  as  the  projection  of  the  current  error  onto  this  space.  In  each 
iteration,  we  carry  out  these  steps  in  the  order  i=  k,  k-1,  ...,  2,  1,  2,  ...  k-1, 
k.  In  [23],  we  establish  the  following  result. 

Theorem  2  The  multiplicative  iterative  refinement  algorithm  has  a  condi- 
tion number  which  is  independent  of  k  and  the  number  of  degrees  of  freedom. 

When  we  now  turn  to  additive  algorithms,  we  note  that  the  most  natu- 
ral algorithm  would  amount  to  using  the  projections  P-,  defined  above,  in 
equation  2.  However,  it  is  relatively  easy  to  show  that  the  condition  num- 
ber of  that  algorithm  grows  linearly  with  k.  A  very  promising  method,  for 
which  to  our  knowledge  the  optimality  has  only  been  established  in  a  quite 
special  model  case  cf.  [14],  is  based  on  using  P-  -  P-^.^  ,  i  <  k  —  1  and  P^  as 
the  basic  projections  in  equation  3.  It  can  be  shown  that  these  differences 
of  projections  are  projections  and  that  the  composite  finite  element  space 
V^  is  the  direct  sum  of  the  corresponding  subspaces.  The  difficulty  experi- 
enced when  trying  to  establish  that  the  sum  of  these  operators  is  uniformly 
well  conditioned  is  related  to  the  lack  of  flexibility  in  representing  a  given 
element  of  F''  as  a  sum  of  elements  of  the  subspaces.  If  we  however  mod- 
ify our  projections,  we  arrive  at  an  algorithm  for  which  we  have  a  proof  of 
optimality. 

Theorem  3  The  additive  method  defined  by  the  projections  P-  -  Pi^2'  *  ^ 
k  —  2,  -P^Jj  and  P^  has  a  condition  number  which  is  independent  of  k  and 
the  number  of  degrees  of  freedom. 

We  note  that  it  is  easy  to  show  that  the  corresponding  subspaces  are 
larger  than  those  of  the  previous  method.  This  makes  our  proof  possible. 
While  the  other  method  is  believed  to  be  faster,  we  have  shown  that  the 
condition  number  of  the  method  of  Theorem  3  is  bounded  by  twice  the 
condition  number  of  the  other.  For  details  see  [23].  We  also  note  that  a 
convergence  proof  was  outlined  in  [5]  for  the  multiplicative  algorithm  in  the 
case  of  k=2  and  that  proofs  and  discussions  are  given  for  both  additive  and 
multiplicative  two  level  algorithms  in  [3]  and  [15]. 
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